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Abstract 



A spin-wave theory of short-range order in the square lattice Heisenberg 
antiferromagnet is formulated. With growing temperature from T = a 
gapless mode is shown to arise simultaneously with opening a gap in the 
conventional spin-wave mode. The spectral intensity is redistributed from 
the latter mode to the former. For low temperatures the theory reproduces 
results of the modified spin-wave theory by M. Takahashi, J. E. Hirsch et al. 
and without fitting parameters gives values of observables in good agreement 
with Monte Carlo results in the temperature range < T < 0.8J where J is 
the exchange constant. 
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Typeset using REVTgX 



Properties of the spin-i quantum Heisenberg antiferromagnet on a square lattice attract 
much attention in connection with the investigation of cuprate-perovskite high-temperature 
superconductors. In accord with the existing theoriesEHi the spectrum of this antiferromag- 
net contains the doubly degenerate (in the magnetic Brillouin zone) magnon mode which 
is gapless at zero temperature. For T > in this mode a gap opens near the center of 
the zone. The appearance of this gap is connected with the short-range antiferromagnetic 
ordering which is established in two-dimensional antiferromagnets at nonzero temperature .0 

Of special note are the spin-wave theory of Refs. |3]|| where without fitting parameters 
values of many observables were calculated in good agreement with the exact diagonalization 
and Monte Carlo calculations in the temperature range < T < 0.6J. Results of these 
works can be obtained with the mean-field decoupling of terms of the Hamiltonian which 
are quartic in the magnon operators or, equivalently, by decoupling of many-particle Green's 
functions in the equations of motion for the one-particle Green's functionsJ3 In the theory 
of Refs. [|,|] the gap in the finite-temperature magnon spectrum appears with imposing the 
constraint of zero staggered magnetization to retain the sublattice symmetry in short-range 
order and to ensure zero site magnetization in the absence of magnetic field. Notice that 
this spin-wave approximation is not rotationally invariant. 

In the present paper we try to obtain a better approximation for the magnon Green's 
functions by transferring the decoupling to higher order equations of motion. This allows 
us to expand the temperature range where the theory conforms with Monte Carlo data up 
to T « 0.8J. Besides, in our theory, as soon as the temperature exceeds zero an additional 
gapless mode arises simultaneously with opening a gap in the conventional magnon mode. 
For large crystals and low temperatures the spectral intensity of the gapless mode is weak 
and our theory reproduces results of Refs. |^,|. With growing temperature the spectral 
intensity is redistributed from the conventional to the gapless mode. 

We consider the antiferromagnetic Heisenberg model on a plane square lattice with the 
Hamiltonian 

H = J^S l S l+a . (1) 

la 

Here Si is the spin-| operator, 1 runs over sites of one of two sublattices, and a are vectors 
of the four nearest neighbors of site zero. In the following discussion the exchange constant 
J is taken as the unit of energy. 

The Dyson-Maleev transformation! is used to represent the spin operators by boson 
operators a\ and b m on the two sublattices A and B 

S{~ = a\, Si~ — (l — a\aij a\, S* = — — a\ai, 
1 G A, 



S m — b m , — (l b^bm 
S z m = ~ + b ] m b m , meB. 

In the new notations Hamiltonian (|l|) acquires the form 



(2) 
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Z 1 m 

+ E ( i a l a l a l & l+a + ^l&l+a & l+a & l+a 



la 



-a\aib\ +a bi +a - -a\b\ +gL - -a x h +i 



(3) 



where N is the total number of sites. 

To investigate the spectrum of elementary excitations we shall calculate the following 
Green's functions: 



£>i(qt) = -i6{t){[a\{t),a^), 
D 2 ( q t) = ~t6(t)([b^(t),a q \), 



(4) 



where a q = (2/iV) 1//2 J2i exp(«ql)e4, fe q = (2/iV) 1 / 2 S m exp(— iqm)b m with the wave vector q 
in the magnetic Brillouin zone, a q (i) = exp(iHt)aJ exp(— iHt), and angular brackets denote 
thermodynamic averaging. Equations of motion for these Green's functions read 



ij t DMt) = -6(f) - 2D 1 (qt) + 2 7q D 2 (qt) + D 3 (qt), 

ij.D%{qb) = 2D 2 (qt) - 2 7qJ D!(qt) + D 4 (qt), 
where 7 q = 4 X)a ex P(2q a ) an d 

D 3 ( q t) = -tf(*)^Ee^([(a[(0O)Wt) - ^(*)ai(*)U*)) ^ 

' la 



(5) 



qm 



( a L+a(*)&L(*)M*) - aL+a(*)am+a(*)M*)) > fl c 



The decoupling of the many-particle Green's functions D3 and -D4 can be carried out at 
this stage. If, like in Ref. 0L an additional constraint of zero staggered magnetization 



5>f -Ess 







£ a l a l + £ 6 m^r 



iV 



(6) 



is imposed by incorporating it with a Lagrange multiplier in Hamiltonian (Q) and if the 
correlations C\ = (ai&i+ a ) = (a\b\ +a ) and (a\a\) = (&m&m) are taken into account in the 
decoupling, results of Refs. |^,| are reproduced. 

In this paper we try to obtain a better approximation. For this purpose we derive 
equations of motion for the functions D3 and D4 and carry out the decoupling in the many- 
particle Green's functions which appear in these equations. If in addition to the mentioned 
correlations we take into account the correlations (ajai +ai +a 2 ) (where ai and a 2 are vectors 
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of nearest neighbors of site zero) and analogous correlations on the second sublattice, we 
find 



i^D a {qt) = 4 (d - I) S(t) + ^D 2 (qt) + 
ijDMt) = 4 (Ci - 7 ,*(t) + f 7 q £>i(q*) 



+ U)( 'i ( ( 'j - -)(l-75) 



A(qt), 



- + 16Ci C a 
2 V 2 



(7) 



£>2(qt), 



where 



K 



8 faaL 



m+ai <2m+a 2 / \ u m+ai u, m+a 2 



)<Qr 



, a 



,>-4C? . 



Owing to the symmetry with respect to translations and rotations the summation in k does 
not depend on m and To derive Eqs. (|7|) we have taken into account that in accord with 
the condition of zero site magnetization in the absence of magnetic field in short-range order 



1 



(a\ai) = (6L & m) = -• 



(8) 



This condition follows also from constraint (||). Substituting Eqs. (0) into Eqs. (|5]) we get 
for the Fourier transforms of Green's functions 



D 2 (qu) 



(AC, - u>) (uj 2 - u 2 2 ) - § [4 (Cx - J) (l - tJ) - U 



4 7q Ci (u 2 - u 2 2 ) - f^7q 



\U) — UJt I [OJ — OJn 



(9) 



where 



OJ 



1,2 



OJ, 



01 



+ oj 02 + « J ± 



01 



(l 



7 C 



01 



^02) 2 + | (^01 - U 02 f + 



OJ, 



02 



16 Cx 



(10) 



7c 



Green's functions (fj|) contain two poles oj\ and oj 2 which correspond to two branches of 
the magnon spectrum. If k is set to zero, the second pole disappears in Green's functions 
and the remaining pole acquires the dispersion of the conventional linear spin-wave theory 



uj\ = 4Ci — 7q) . For K^fla gap of the width k 1 ^ 2 opens in this branch near q = 0. 
Simultaneously the second branch acquires finite spectral intensity which is subtracted from 
the intensity of the first branch. As will be seen below, in the considered temperature 
range C\ is close to ~. This allows one to approximate the dispersion of the first branch 



as o>i 



K + 16C 2 (l - 7 2 



1/2 



which, with some change of notations, coincides with the 
dispersion obtained in Refs. |||]. Technically in those works the gap appears in the magnon 
spectrum with imposing the constraints of zero staggered or site magnetization, which are 
incorporated into the Hamiltonian or free energy. In the more exact treatment of the Green's 
functions the magnon gap arises without such changes of the Hamiltonian. 
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Parameters G\ and k in the above equations contain correlations which can be deduced 
from Green's functions (||). This gives the following self-consistency conditions for evaluating 
the parameters: 



q 



(11) 



where K x = (2/N) E^cos(q x - q y ), K 2 = (2/N) £ q T q cos(2 ga: ) 



(12) 



J* 



uox — U), 



02 



J, 2« (d - ^) (l 



(13) 



7^ 



Tq 



1 



n(a>i) + 



1 

CJ 2 



n(o; 2 ) + 



n(c<j) = [exp(u/T) — 1]~ , and q y are the components of the wave vector q (the intersite 
distance is taken as the unit of length). In Eqs. (|TTD, (|T3) and below summations over wave 
vectors are carried out over the magnetic Brillouin zone. Analogously, from Eqs. (H) we find 
for the spin correlation functions 



(SiS r ) = 
(SiS m ) 



1 2 



AT 2 - q 
q 



AT 



D iq(l-m) 



(14) 



It can be seen that correlations (£i + <Sw m )) are zero and only (S{Sy^) contribute in the above 
expressions. Thus the considered spin- wave approximation is not rotationally invariant. 

To determine the parameters C\ and k, instead of one of Eqs. (Pj), (0) one can use the 
condition of zero site magnetization (Rl) which can be rewritten as 



N ^ q 

q 



(15) 



Comparing results obtained with different pairs of these three equations with Monte Carlo 
data, we found that the best agreement is achieved with Eqs. ([□]) and ([TBI). Notice also 
that all three equations can be used introducing some parameter in addition to k and 
C\. In particular, such additional parameter appears if we incorporate the constraint of 
zero staggered magnetization (|6|) in Hamiltonian (|3[) with a Lagrange multiplier [condition 
(|15|) follows from this constraint, if one takes into account the translation and sublattice 
symmetry] . This approach is a generalization of the method of Refs. |3],|] — results of these 
works can be obtained from formulas of the three-parameter approach by setting k = 0. We 
postpone the consideration of this approach to the end of the article. For now we discuss 
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results obtained in the two-parameter approach based on Eqs. (|TT|) and (|i~5|). Notice that 
the relation (S 2 ) = f follows from Eqs. fll4]) and C|T5|) . 

Let us first consider the case of low temperatures. As will be seen below, k = 0(N~ 2 ) 
for T = 0. Since q 2 = (2ir) 2 (n 2 + n 2 )/N where n x and n y are integers, even for the smallest 
q 2 , excluding q 2 = 0, and Uq 2 are much larger than k for large N. In this case Eqs. ( |T3"D 
can be simplified to 

>4 



4d 






n 






4C l7q 





n{u x ) + - 



(16) 



O(zc). 



where uo\ 



[OJ, 



01 



+ ft) 1 / 2 (here we took into account that C\ — | <C Ci). Notice that the 



contribution of the gapless magnon mode dropped out from these equations. They are also 
suitable for an infinite crystal when T — > and k is exponentially small. With these Z q , J7" q 
and cui Eqs. (|Tl"l ) and ([15]) come close to the respective formulas of Ref. [| As a consequence, 
in the limit iV — > oo, T — > values of observables obtained with Eqs. (|ll]), (HH) are similar 
to those found in Ref. |3[. Therefore we only briefly discuss this limit below. 
Using Eqs. (|16|) , for large N and T 

as 



iV 



3 jqr 



4Ci 



the sum in the first equation ([14]) can be written 

(17) 



+ 



1 



N 



E 



i 



7 2 

'q 



The sublattice magnetization m 



AT 



1 Eq^O (l 



7 



-1/2 



4Ci/(Nk 1 / 2 ) is determined by Eq. 



For an infinite crystal m 



0), m = 1 - 
0.3034 which is in good agreement 



with the Monte Carlo calculations.l-H--l As follows from the above formulas, k is actually of 
the order of N~ 2 . From Eq. (|T?D for large r we get m + (2 l / 2 nr^ and after analogous 



transformations in the second equation ([14]) we find 

\/27rr 



(S S r ) 



m + 



-i 



(18) 



where (— l) r = +1 or —1 depending on whether the sites and r belong to the same or 
different sublattices. In Table |the zero-temperature spin correlations C r = (S S r ), obtained 
by numerical solution of Eqs. (O) and ( 15 ) (see below), are compared with results of the 
projected Monte Carlo method.^ The values agree nicely. 
For iV — > oo, T — > we find from Eqs. ([TT|), ( |T5D and 

C\ = mi 



4 C(3) (^-) +0(T S 



7T 



Texp 

2 \l/2 



\4mi 
f 



(19) 



1 + 0{T 2 



where m\ = 1 — N 1 ^ q y- — 7 q J ' = 0.57897 and C{ x ) is the Riemann zeta function. For 
large r the sums in the spin correlations (SoS r ) flUP can be rewritten as 

T ff d 2 q ,„_ T 



3 «qr 



2vr) 2 C 1 77 g 2 + (20 " 2 
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with the correlation length 

£ = Cl f- = ^ exp [j + op-)] . (20 ) 

This value of the correlation length is in agreement with results obtained in Refs. [l],|2|. 

To solve Eqs. QTT]) and (|T5|) for arbitrary T and N we determined the minimum of the 
function 

^„ K ) = (i-|x:i q ) 2 + (i-^e^) 2 , 

which is constructed from squares of the differences of the right and left sides of these 
equations. The iteration procedure with the steepest descent method was continued until 
F was less than 10~ 12 . As an example, parameters obtained by this procedure for a 20x20 
lattice are given in Table [Tl| These parameters were used for calculating the static uniform 
susceptibility 



i 

3T 



N ^ ^ J ^ 4 



(21) 



and the energy per spin 



E = 2(SiSi +a ) = -2Cl (22) 



Results for a 20x20 lattice are shown in Figs. |l] and ||] together with the data obtained in 
the Monte Carlo calculations0@ and in the modified spin-wave theory of Refs. 0,[l[ As seen 
from the figures, results obtained in the spin-wave approximations of the present work and of 
Refs. ||,f| are close and are in good agreement with the Monte Carlo results in the considered 
temperature range. The size dependence of \ and E calculated with the above formulas is 
negligible for large enough lattices — the difference between values obtained for a 40x40 
lattice and those shown in Figs, [I] and is less than the size of symbols in these figures. 
Notice also that zero frequency of the second branch at q = does not lead to divergencies 
in the above formulas as the respective numerators in Eqs. (|1^) approach zero sufficiently 
rapidly when q — > 0. 

As seen from Table [II], starting from low temperatures the gap determined by k grows 
with T. With this growth the spectral intensity of the second branch and its contribution 
to the observables increases. However, near T = 0.5 k starts to decrease. For T « 0.7, when 
k — ► and C\ ~ 0.5, the considered two-parameter approximation ceases to work — it is 
impossible to find parameters for which Eqs. flTTp and (|T^) are fulfilled with high accuracy. 

Now let us consider the three-parameter approximation mentioned above. The additional 
parameter is the Lagrange multiplier A with which the terms J2i a i°i + Em ^L^m are added 
to Hamiltonian (0) to ensure the fulfillment of constraint (|6]). With these additional terms 
Di(qu) in Eq. (^), k in Eq. (|T2|) and X q in Eq. flT3|) are slightly changed 

= (^y- 1 - u) (c 2 - c 2 2 ) - | [4 (Cx - \) (l - t») - u 



7 



1 - 



?77q 



7 2 
'q 



1 



7l(Wi) 



1 

CJ 2 



n(o> 2 ) + 



whereas /^(qu;) in Eq. (| 
In these equations 77 = [1 



Jq in Eq. 

A/(4C 1 )]- 1 



, Eqs. (|n|), (|T3D, (gT]) and (0) retain their form. 



1,2 



01 



+ ^02 + K J 



±4w 



01 



i6c? 



1 



02) 



7 2 
'q 



+ 



ft 



01 



02 J 



02 



16d (C; 



1 



7; f (1 



To solve Eqs. flTTD, (p.2f) and ( 15 ) we determined the minimum of the function F'(Ci, ft, 77), 
constructed from squares of the differences of left and right sides of these equations, by 
the steepest descent method. The static uniform susceptibility and the energy per spin 
calculated with the parameters obtained in this way for a 20x20 lattice are shown in Figs. [I] 
and |^. As seen from the figures, in comparison with the two-parameter approximation 
the three-parameter approach agrees slightly better with the Monte Carlo results and is 
applicable in the wider temperature range < T < 0.8. However, the convergence of the 
steepest descent method in the three-parameter approximation is much worse than in the 
two-parameter one. We connect this with the observation that both approaches give similar 
pictures of magnon modes and in the three-parameter approximation the two parameters k 
and A (or rj) determine one physical quantity — the magnon gap. In this approximation the 
worsened convergence is connected with the flat minimum of F' considered as the function 
of k and A. 

In the present work we have found the gapless mode in the state with short-range order. 
In ordered states excitations of such type correspond to Goldstone's mode and point to the 
existence of the continuous degeneracy of these states.^ The observation of the gapless mode 
in the considered disordered state may also be connected with the continuous degeneracy of 
this state. In some respects a similar magnon spectrum of short-range order was obtained in 
Refs. [iiyr5| . When considered in the magnetic Brillouin zone the spectrum consists also of 
two modes one of which has a gap near q = and another is gapless. However, the shape of 
these branches, their relative spectral intensities and temperature behavior differ essentially 
from those obtained here. The spin- wave approximations of Refs. H,15 are rotationally 
invariant but are not self-consistent — correction coefficients and data of other calculations 
are needed to fit results to experiment in the temperature range T < 1. 

In conclusion, the spin-wave theory of short-range order in the square lattice Heisenberg 
antiferromagnet was formulated. In agreement with previously obtained results we found 
that as soon as the temperature exceeds zero and long-range order gives way to short-range 
order a gap opens in the conventional magnon mode near q = 0. We found additionally 
that a new gapless mode arises simultaneously with opening the gap. With growing tem- 
perature the spectral intensity is redistributed from the conventional to the gapless mode. 
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We considered spin-wave approximations with two and three self-consistently determined 
parameters. For low temperatures and large crystals the theory reproduces results of the 
modified spin-wave theory of Refs. |3|,|]. Without fitting parameters our calculations give val- 
ues of the static uniform susceptibility and energy per spin in good agreement with Monte 
Carlo results in the temperature range < T < 0.8J where J is the exchange constant. 
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FIGURES 



FIG. 1. The static uniform susceptibility obtained in the Monte Carlo simulation for a 12x12 
latticed (•), in the modified spin- wave approximation of Refs. |3],|| (o) and in the two-parameter 
(+) and three-parameter (x) spin- wave approximations of the present work. In the spin- wave 
calculations a 20x20 lattice was used. 

FIG. 2. The energy per spin obtained in the Monte Carlo simulationEl (•), in the modified 
spin- wave approximations of Refs. ||,|| (o) and in the two-parameter (+) and three-parameter (x) 
spin- wave approximations of the present work. In the spin- wave calculations a 20x20 lattice was 
used. 
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TABLES 



TABLE I. The zero-temperature spin correlations C r obtained with the two-parameter 
spin- wave approximation (SW) for a 20x20 lattice in comparison with the projected Monte Carlo 
data (PMC)fl 





PMC 


SW 


Ci,o 


-0.3348 


-0.3354 




0.2028 


0.2016 


C2,0 


0.1772 


0.1751 


C2,l 


-0.1671 


-0.1648 


C2,2 


0.1475 


0.1454 


C3,0 


-0.1491 


-0.1461 


C3,l 


0.1430 


0.1404 



TABLE II. Parameters C\ and k obtained from Eqs. flTl| ) and ([[!]) for a 20x20 lattice. 



T 


Ci 


K 


10" 5 


0.57912 


2.773 • 10" 4 


0.1 


0.57904 


4.685 • 10" 4 


0.2 


0.57826 


9.390 • 10~ 4 


0.3 


0.57543 


1.443 • 10~ 3 


0.4 


0.56819 


1.941 • 10~ 3 


0.5 


0.55409 


2.303 • 10~ 3 


0.6 


0.53037 


2.093 • 10~ 3 


0.7 


0.50003 


3.664 • 10" 6 
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Fig. 1. A. Sherman and M. Schreiber, Temperature behavior 




2. A. Sherman and M. Schreiber, Temperature behavior 



